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Abstract 

In this paper, the connections are investigated between two different approaches towards 
the parametrization of multivariable stable all-pass systems in discrete-time. The first ap- 
proach involves the tangential Schur algorithm, which employs linear fractional transforma- 
tions. It stems from the theory of reproducing kernel Hilbert spaces and enables the direct 
construction of overlapping local parametrizations using Schur parameters and interpolation 
points. The second approach proceeds in terms of state-space realizations. In the scalar case, 
a balanced canonical form exists that can also be parametrized by Schur parameters. This 
canonical form can be constructed recursively, using unitary matrix operations. Here, this 
procedure is generalized to the multivariable case by establishing the connections with the 
first approach. It gives rise to balanced realizations and overlapping canonical forms directly 
in terms of the parameters used in the tangential Schur algorithm. 



1 Introduction 



Stable all-pass systems of finite order have several applications in linear systems theory. Within 
the fields of system identification, approximation and model reduction, they have been used 
in connection with the Douglas-Shapiro-Shields factorization, see e.g., [U [TO1 HH] , to obtain 
effective algorithms for various purposes, such as H2— model order reduction of stable linear 
systems. The class of stable all-pass transfer functions of finite order is bijectively related to 
the class of rational inner functions, of which the differential structure has been studied in pQ. 
There, a parametrization for the multivariable case has been obtained by means of a recursive 
procedure, the tangential Schur algorithm, which involves Schur parameter vectors, interpolation 
points and normalized direction vectors. In the scalar case, a single coordinate chart suffices to 
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entirely describe the manifold of stable all-pass (or lossless) systems of a fixed finite order. In 
the multivariable case, the approach leads to atlases of generic charts covering these manifolds. 

In another line of research, balanced state-space canonical forms have been constructed for 
various classes of linear systems, with special properties of these classes (such as stability) built 
in, see e.g., |2H 116]. Balanced realizations are well-known to have numerical advantages and are 
useful for model reduction purposes in conjunction with balance-and-truncate type procedures. 
A balanced canonical form for SISO stable all-pass systems in continuous time was presented 
in [22]. In the constructions of |15| . [16] . the case of stable all-pass systems in continuous- 
time plays a central role. In the scalar case, the resulting canonical form for lossless systems 
is balanced with a positive upper triangular reachability matrix. In the multivariable case, 
Kronecker indices and nice selections are used to arrive at balanced overlapping canonical forms 
for lossless systems. For discrete-time stable all-pass systems, canonical forms can be obtained 
from the results in continuous-time by application of a bilinear transformation. However, this 
destroys certain nice properties of the canonical form; e.g., truncation of state components no 
longer leads to reduced order systems that are balanced and in canonical form. Therefore, the 
ideas of |15| and [16] are applied in |14] to the scalar discrete-time stable all-pass case directly. 
This leads to a balanced canonical form with the desired properties, for which it turns out that 
it can in fact be parametrized using Schur parameters. 

In this paper, the connections between these two approaches are investigated and the results of 
[2] are generalized to the multivariable case. Moreover, the atlases developed in [1] are supplied 
with balanced state-space realizations in terms of the parameters used in the tangential Schur 
algorithm. In Section [31 we introduce a mapping which acts on proper rational matrix functions 
and which has a particularly simple expression on the level of realizations, see formula ()34p . 
When specialized to the case of lossless functions, it allows for a recursive construction of bal- 
anced realizations that can be implemented as a product of unitary matrices. This construction 
has several advantages from a numerical point of view. In Section HI we recall some basic results 
on linear fractional transformations of lossless functions, which are at the heart of the tangential 
Schur algorithm, described in Section [5j The results of these sections are mostly well-known, 
except for the nice factorization of J-inner functions in Proposition 15. 2\ which we could not find 
in the literature and which enables an elegant presentation of these results. The main technical 
results are presented in Section [6] and they describe the complete connection between the con- 
struction of Section [3] and the tangential Schur algorithm in its more general form. In Section [7] 
the results are used to construct atlases for the manifold of lossless systems of a fixed McMillan 
degree and for the manifold of stable linear systems of a fixed McMillan degree. 



2 Preliminaries 



In this section a brief overview of results from the literature is presented to provide the necessary 
background for the main constructions of this paper. They concern realization theory, balancing 
and the basic theory of lossless and J-inner functions. 
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Realization theory and balancing. 

Consider a linear time-invariant state-space system in discrete time with m inputs and p outputs: 



x t +i = Ax t + Bu t , 
yt = Cxt + But, 



with t G Z, xt G C n for some nonnegative integer n, Ut G C m , G C p . Furthermore, the 
matrices A, B, C and B with complex-valued entries are of compatible sizes: n x n, n x m, 
p x n and p x m, respectively. The corresponding transfer matrix of this system is given by 
G(z) = B + C{zl n — A)~ 1 B, which isapxm matrix with rational functions as its entries. Two 
state-space systems are called input-output equivalent if they have the same transfer matrix. 
Here, two rational (matrix) functions are identified if they agree almost everywhere, i.e., common 
factors are always canceled. Note that a transfer matrix in this set-up is always proper (i.e., 
analytic at infinity), B being the value at infinity. 

From realization theory it follows that, conversely, any p x m rational matrix function G(z) 
analytic at infinity can be written in the form 



G{z) = B + C(zI n - A)- l B 



(3) 



where (^4, B, C, B) is an appropriate quadruple of matrices and n a suitable state space dimen- 
sion. Such a quadruple with the associated expression ([3]) is called a state-space realization of 
G{z). To such a realization we associate the block-partitioned matrix 



R 



B C 
B A 



(4) 



which we call the realization matrix as in J14J. It will play an important role in the sequel. 



If all possible realizations of G{z) have state-space dimension at least as big as the state-space 
dimension n of (A, B, C, B) then the latter realization is said to be minimal, and n is called the 
order or the McMillan degree of the transfer function. Two minimal realizations (A, B, C, B) 
and (A' , B' ,C' , B') of a given function G(z) are always similar: there exists a unique invertible 
matrix T such that 



(5) 



" B' 


C' ' 






" 




' B 


C ' 




'Ira 


B' 


A' 







T 




B 


A 




T- 1 



As is well-known, an output pair (C, ^4) is observable if the observability matrix 

C 

CA 



(6) 



CA™ 1 ' 1 

has full column rank. An input pair (A, B) is reachable if the associated reachability matrix 

[B AB ... A n ~ 1 B } (7) 
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has full row rank. It is well-known that a realization (A, B,C,D) is minimal if and only if (C, ^4) 
is observable and (A, B) is reachable; cf., e.g., [IT]- In this case, the poles of the function G(z) 
are the eigenvalues of A, and its McMillan degree is equal to the sum of the degrees (see [EJ Sect. 
6.5]) of all the poles of G(z). This provides an alternative definition for the McMillan degree, 
which generalizes to the larger class of proper and non-proper rational functions by including a 
possible pole at infinity with its appropriate degree. 

Let (A,B,C,D) be some realization of a transfer function. If the eigenvalues of A all belong 
to the open unit disk, then the matrix A is called (discrete-time) asymptotically stable, and 
(A, B, C, D) an asymptotically stable realization. If (A, B, C, D) is an asymptotically stable 
realization, then the controllability Gramian W c and the observability Gramian W are well 
defined as the exponentially convergent series 

oo oo 

W c = Y J AkBB *( A *) k , W = Y J ( A *) k C*CA k , (8) 

k=0 k=0 

where the notation * is used to denote Hermitian transposition of a matrix (i.e., the joint action 
of matrix transposition and complex conjugation of the matrix entries). The Gramians are 
characterized as the unique (and positive semi-definite) solutions of the respective Lyapunov- 
Stein equations 

W C -AW C A* = BB*, (9) 
W -A*W A = C*C. (10) 

Moreover, under asymptotic stability of A it holds that W c is positive definite if and only if the 
pair (A, B) is reachable, and W Q is positive definite if and only if the pair (C, A) is observable. 
A minimal and asymptotically stable realization (A, B, C, D) of a transfer function is called 
balanced if its observability and controllability Gramians are both diagonal and equal. Any 
minimal and asymptotically stable realization (A, B, C, D) is similar to a balanced realization. 
The concept of balanced realizations was first introduced in [20J in the continuous time case 
and used for model reduction. In [25] the same was done for the discrete time case. Balanced 
realizations are now a well-established tool which often exhibit good numerical properties. 

J -inner, J -unitary and lossless functions. 

For any matrix function R(z), we define the matrix functions R*(z) and R)(z) by 

R*(z) := R{z)\ and R\z) := R*( Z - 1 ). (11) 

Note that if z lies on the unit circle, then W(z) = R{z)* . For (square) Hermitian matrices P 
and Q, either of the notations P < Q and Q > P will be used to express that Q — P is a positive 
semi-definite matrix. As is well-known this introduces a partial ordering on the set of square 
Hermitian matrices. 



Now let S be a constant k x k matrix which is both Hermitian and unitary. Note that E is 
unitarily similar to a signature matrix, i.e., there exists a unitary k x k matrix U for which 



UT,U* attains the form 









-It 



for some non- negative integers q and r with q + r = k. 
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A square rational matrix function Q(z) of size k x k is called E-inner (in the unit disk), if at 
every point of analyticity z of it satisfies 



e(z)*se(z) < e, |z| < 1, 
e(z)*se(z) = s, |z| = i, 
e(z)*se(z) > s, |z| > i. 



(12) 
(13) 
(14) 



It can be shown that this definition contains redundancy in the sense that either one of these 
three defining properties (fT2l) - (fT4l) is implied by the other two; cf., e.g., [23| [8l [T2]. 

Any rational matrix function @{z) which satisfies the property (|13p is called E-unitary. For such 
functions, replacing z by z^ 1 on the unit circle, the identity (I13p extends almost everywhere by 
analytic continuation, so that a rational E-unitary function is invertible and its inverse is given 



This applies a fortiori to E-inner functions. A function is called E-lossless if it is (— E)-inner. If 
Q(z) is E-inner, then Q*(z) is also E-inner whereas both ©"(#) and 0(z) -1 are E-lossless. The 
class of E-inner functions is closed under multiplication. 

As usual, an /^-lossless function is just called (discrete-time) lossless, or (discrete-time) stable 
all-pass and an I^-inner function is simply called inner. Throughout this paper we will be 
much concerned with lossless and inner functions of size p x p, and with J-unitary and J-inner 
functions of size 2p x 2p, where J denotes the following signature matrix: 



Note that J-inner and J-lossless functions in general may have poles everywhere in the complex 
plane, but inner functions are analytic inside the unit disk and lossless functions are analytic 
outside the unit disk (including at the point at infinity). Therefore a rational lossless function 
G(z) is proper, having an inverse which is inner and given by 



A scalar rational function is a Blaschke product if it maps the unit circle into the unit circle. It 
can be shown that such a function can be written as the product of a finite number of Blaschke 
factors, which are rational functions of first order that map the unit circle into the unit circle. 

Lemma 2.1 Let G(z) be a p x p lossless function of McMillan degree n. Then its determinant 
detG(z) is a Blaschke product. If one writes it in the form det G(z) = e{z)/d{z) with e and d 
co-prime polynomials, then the degree of d(z) is equal to n. 

Proof. See e.g. [5]. □ 



by 




(15) 




(17) 
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We finally describe some particular J-inner functions that will be intensively used in the sequel, 
namely the constant J-unitary matrices and the J-inner functions of McMillan degree one. 
Details and proofs can be found in [8]. 

Every constant J-unitary matrix can be represented in a unique way (see [H Thm. 1.2]) as 
follows: 

" P 



M = H{E) 



Q 



(18) 



where P and Q are p x p unitary matrices and H{E) denotes the Halmos extension of a strictly 
contractive pxp matrix E (i.e., such that I — E*E > 0). This Halmos extension H(E) is defined 
by 



H(E) 



(I -EE*)' 1 / 2 E(I - E*E)- 1 / 2 
E*(I - EE*)- 1 ' 2 (1-E*E)- 1 I 2 



(I -EE*)- 1 / 2 (J — EE*)~ X / 2 E 
(I — E* E)~ 1 / 2 E* (I-E*E)- 1 / 2 



(19) 

It holds that H(E) is Hermitian, J-unitary and invertible with inverse H{E) 1 = H(—E). Also, 



KH(E)K = H(E*) for K 



I p 
I P 



Two other important forms of J-unitary matrices are described in the following lemma ( see [HI 
Thm.1.3]). 



(22) 



Lemma 2.2 If x is a 2p-vector such that x*Jx ^ 0, then the set of2p x 2p matrices of the form 

X x (a) = I 2p + {a-l)x(x*Jx)- 1 x*J, a G C, (20) 

is closed under multiplication. We have that 

X x (a)X x ((3) = X x (af3), detX x (a) = a, X x {a)~ l = X x {a~ x ), for a ^ 0, (21) 

J - X x (a)JX x ((3)* = (1 - a(3*)x(x*Jx)~ 1 x*J, 
and X x (a) is J-unitary if and only if \a\ = 1. 

If x is a 2p-vector such that x*Jx = 0, then the set of matrices of the form 

Y x (a) = I2 P + axx*J, a £ C, 

is closed under multiplication. We have that 

Y x {a)Y x {P) = Y x (a + (3), detY x (a) = 1, Y^a)' 1 = Y x (-a), 

J - Y x (a)JY x ((3)* = (a + (3*)xx*, 
and Y x (a) is J-unitary if and only if a is pure imaginary. 



(23) 

(24) 
(25) 
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We now come to the description of J-inner matrix functions of McMillan degree 1. We shall call 
these functions 'elementary J-inner functions' or 'elementary J-inner factors', since it has been 
shown by Potapov (see [23]) that every J-inner matrix function of McMillan degree n can be 
decomposed into a product of n such elementary J-inner factors. Elementary J-inner functions 
are obtained by applying the functions X x and Y x to appropriate scalar functions of degree one, 
namely the Blaschke factor b w {z) defined by 



and the Caratheodory function, 



where T denotes the unit circle. 



b w {z):=^-^, wiT (26) 
1 — wz 



7 A- 11) 

cu,(z):=^—, weT, (27) 
z — w 



Theorem 2.3 Let $(z) be an elementary J-inner factor with a pole at z = 1/w. Then, apart 
from a constant J -unitary multiplier on the right, &(z) must be in one of the following three 
forms: 

(1) Ifw G C, w $ T, 

<&(*) = X x (b w {z)) = + (b w (z) - 1)^, (28) 
(see Eqn. i20\) ) for some x G C 2p such that : 

7^T12 > °- ( 29 ) 
1 — \wr 



(2) IfweT, 

= Y x (-5 c w {z)) =h p -6 c w (z) xx*, J, (30) 
(see Eqn. $23\)) for some nonzero x G C 2p such that x* Jx = 0, and for some real 5 > 0. 

(3) If w = co, then 

. W =*(i)-/» + (I-l)pg, (31) 



for some x G C p such that x* Jx < 0. 



3 Recursive construction of balanced realizations of lossless sys- 
tems 



In this section we will present a number of results that will be used in the proposed recursive 
construction of balanced state-space realizations of lossless functions. 
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With each pair (U, V) of (p + 1) x (p+ 1) matrices we associate a mapping Tjjy which is defined 
to act on proper rational p x p matrix functions G(z) as follows. 



T uy : G{z) ^ F x (z) + 



F 2 (z)F 3 (z) 



z-F 4 (z) ' 

with F\{z) of size p x p, F 2 {z) of size p x 1, F 3 (z) of size 1 x p and F±(z) scalar, defined by: 



(32) 



F(z) 



F x {z) F 2 (z) 
F 3 (z) F 4 (z) 



V 



1 

G(z) 



U*. 



(33) 



A state-space realization of G(z) = J r u,v(G(z)) can be obtained by working directly on any 
realization matrix R of G(z). This is the content of the following proposition. 



Proposition 3.1 Let G{z) be a p x p proper rational transfer function and (U,V) a pair of 
(j) + 1) x (p + 1) matrices, then G(z) = Tu,v(G(z)) is well-defined. Let (A,B,C,D) be a 
state-space realization of G(z) with n- dimensional state-space. Then a state-space realization 
(A, B,C , D) of G(z) with (n + 1)- dimensional state-space is given by: 



D 


C 




' V 





B 


A 







In _ 



1 
D C 
B A 



U* 

In 



(34) 



Proof. Since G(z) is proper, F^{z) is proper too and therefore z — F±{z) does not vanish identi- 
cally so that Tjjy{G{zj) is well-defined. 



Now observe that the right-hand side of Eqn. (I34p provides a realization matrix with n-dimensional 
state-space for the transfer function F(z) given by Eqn. (j33|) : 





C F ' 




' V 


B F 


A F 




I n 



1 
D C 
B A 



U* 








Here Dp is {jp + 1) x (p + 1), Cp is (p + 1) x n, Bp is n x (p + 1) and = A is n x n. We shall 
prove that changing the partitioning of this matrix to the one in ^54 



' Dp 


C F ' 




D 


C 


Bp 


A F 




B 


A 



where D is p x p, C is p x (n + 1), B is (n + 1) x p and A is (n + 1) X (n + 1), gives a realization of 
Gr(z) with n + 1-dimensional state-space. To see this, let the matrices C, B and ^4 be partitioned 
as follows: 



D C 
B A 



D 


Ci 


c 2 




D 


Ci 


c 2 


#1 








Bi 




A 2 


B 2 


^3 


Af 




B 2 


A 3 


A_F 
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where C\ is p x 1, C 2 is p x ra, Si is 1 X p, £? 2 is n x A 2 is 1 x ra, ^3 is n x 1 and A 4 is scalar. 
Note that it follows that F 4 (z) has the realization (Ap, A3, A 2 , A±) and therefore 

z-Ax- A 2 {zl n - A F )~ 1 A 3 = z- F 4 (z). 

Using this together with the well-known formula for the inverse of a partitioned matrix (see [8] 
or [T7]) we compute 



(z/, 



n+l 



4)" 



(z - F4JZ))- 1 -(z - ^(z))" 1 ^^ - Af)" 1 

-(zln-A^Asiz-F^z))- 1 A(z) 



with 



A(z) = (z/„, - A F )~ 1 + (z/ n - ^r%(^ - F 4 (z))- 1 A 2 (z/„ - A F )-\ 
We then have 



D + C(zl n+ i — A)~ l B 



D + 



Ci Co 



(zl, 



n+l 



AT 1 



By 
B 2 



= D + Cxiz- F^z))- 1 ^ + C x (z - F 4 (z))- 1 2 2 (zJ n - A F y l B 2 + 
+C 2 (z/ n - Af)- 1 ^^ - F±{z))- l B\ + C 2 (zl n - Af)' 1 ^ + 
+C 2 (zI n - A F y l A 3 (z - F 4 (z))- l A 2 (zI n - A F y 1 B 2 

= F 1 (z)+F 2 (z)(z-F 4 (z))- 1 F 3 (z), 

which finally proves that (A, B, C, D) is a state-space realization of G{z) with (n+l)-dimensional 
state-space. □ 



The following proposition characterizes the nature of minimal balanced state-space realizations, 
in discrete time, of rational lossless functions. The fact that a lossless function admits a minimal 
balanced realization for which the realization matrix is unitary is well-known from the literature; 
see, e.g., [12J. However, the converse result, which asserts that unitary realization matrices 
correspond to possibly non-minimal realizations of lossless functions, appears to be novel and 
generalizes |14| Prop. 2.4] to the multivariable case. 



Proposition 3.2 (i) For any minimal balanced realization of apxp rational lossless function the 
observability and controllability Gramians are both equal to the identity matrix and the associated 
realization matrix is unitary. 

(ii) Conversely, if the realization matrix associated with a realization (A,B,C,D) of order n 
of some p x p rational function G is unitary, then G is lossless of McMillan degree < n. The 
realization is minimal if and only if A is asymptotically stable and then it is balanced. 



Proof, (i) As is well known, any minimal realization (A, B,C, D) of a lossless function G is 
asymptotically stable and satisfies (see [T2] ) 



" D 


C ' 




' h 


" 




- D * 


B* ' 






" 


B 


A 







p 




c* 


A* 







p 
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for some unique positive definite n x n matrix P which is precisely the controllability Gramian. 
Since P is positive definite, it can be reduced to identity by using any state space trans- 
formation matrix T such that P = TT*, producing a new minimal state space realization 
(A,B,C,D) = {T- X AT,T- X B,CT,D). Then Eqn. ([35} asserts that the associated new real- 
ization matrix is unitary. Upon changing the order of multiplication of this unitary realization 
matrix and its Hermitian transpose it follows that the new observability Gramian is also equal 
to identity. Hence, ( A, B, C, D) is balanced. Any other minimal balanced realization of G{z) 
can be obtained from (A, B, C, D) by means of a unitary change of basis of the state space, 
which leaves the associated realization matrix unitary. 

(ii) Conversely, if a realization (A, B, C, D) has unitary realization matrix, a standard straight- 
forward calculation shows that 

I p - G(z)G{w)* = (zw - l)C(zI n - A)-\wI n - A*)- X C*, (36) 

for every choice of z and w outside the spectrum of A. Setting w = z, it follows that G is lossless. 
Therefore, G has all its poles inside the open unit disk, which makes asymptotic stability of A 
a necessary condition for minimality of (A, B,C, D). On the other hand asymptotic stability 
of A implies that minimality of (A, B,C, D) is equivalent to the existence of positive definite 
solutions to the Lyapunov-Stein equations ([9])- (|10p . Unitarity of the realization matrix, however, 
implies that I n is always a positive definite solution to both these equations, showing asymptotic 
stability of A to be a sufficient condition for minimality of (A, B,C,D) as well. Clearly in such 
a case (A, B ,C, D) is also balanced. □ 



Remark. Let G\ and G2 be two lossless functions of McMillan degree n\ and ri2, respectively, 
with minimal realizations (Ax, Bx, Cx, D\) and (A2, B2, C2, D2) both having unitary associated 
realization matrices. Then the cascade realization [A, B, C, D) of the lossless function G = G1G2 
is obtained as 



" D 


C ' 


B 


A 








Ci " 










Bi 





A x 



D 2 


c 2 







\ DxD 2 


DxC 2 


Cx 


B 2 


A 2 







B 2 


A 2 















BxD 2 


-B1C2 


Ax 



which yields again a unitary realization matrix which is minimal. 



Three properties of the class of mappings J-uy are collected in the following lemma for later 
use. These properties follow straightforwardly from the definition of the mapping Tjjy with U, 
V arbitrary (p + 1) x (p + 1) matrices. 



Lemma 3.3 Let U and V be arbitrary (p + 1) x {p + 1) matrices and P and Q both p x p 
matrices. Then we have: 



(i) T 



1 
P 



y 



1 

Q 



(G(z)) = T V y{QG{z)P*), 



(ii) G(z) = F uy {G{z)) ^ G*(z) = F VJU {G*{z)), 
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(iii) T 





" 


u, 




' 





t _ 





Z _ 



(G(z)) = F uy {G{z)), for all £ E C with |f | = 1. 



V 



The next property shows some of the importance of the J~u,v mappings in relation to lossless 
systems: if U and V are unitary matrices then lossless systems are mapped into lossless systems. 
The McMillan degree increases by at most one under this mapping. 



Proposition 3.4 Let U and V be unitary (p + 1) x (p-f 1) matrices. 

The mapping J-jjv sends a lossless function of order n to a lossless function of order < n + 1. 

Proof. Since G{z) is lossless of order n it has a minimal realization with a unitary realization 
matrix of size n + p. Formula (|34p gives a realization of G = J~ijy{G) which is unitary of size 
n + p + 1. Therefore, G(z) is a lossless function of order at most n + 1. □ 

Starting with a p x p unitary matrix Go, interpreted as a lossless transfer matrix of order 
zero, and applying the previous proposition repeatedly, one obtains a class of lossless transfer 
matrices and corresponding realizations. The class that is obtained after n such recursion steps 
is parametrized by the finite sequence of unitary matrices Ux, V{, U2, V2, ■ ■ ■ , U n , V n , together 
with the unitary matrix Gq. It consists of lossless systems of order < n + 1. In Section [6] the 
question which lossless transfer matrix functions can be obtained in this way will be treated. 
This generalizes the results of [2] for the SISO case, where Uk = h and is a 2 x 2 Householder 
reflection matrix, for k = 1, . . . , n. 



4 Linear fractional transformations 



Linear fractional transformations occur extensively in representation formulas for the solution 
of various interpolation problems [3]. They are at the heart of the parametrization of lossless 
functions through the tangential Schur algorithm. The properties of a linear fractional trans- 
formation of matrices were studied in a form adapted to the needs of J-theory in [24] and in an 
algebraic setting in [26] . For more details on and certain proofs of the basic results presented in 
this section, we refer to these papers. 

Let E C 2px2p (z) be an invertible rational matrix in the variable z, block partitioned as 

"01 02 " 







e 3 0, 



with blocks 0» (i = 1, . . . , 4) of size p x p. 



Associated with 0, let the linear fractional transformations 7e and 7e be defined to act on 
rational matrices G E C pxp (z) in the following way: 

r©:G^(0 4 G + 3 )(02G + 1 )- 1 , (37) 
fe:G^(G02 + 04) _1 (G0 1 + 3 ). (38) 
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The domains of these mappings are denoted by A4q an d -Mo, respectively. They consist of those 
G for which the expressions Q2G + ©i and G@2 + @4 ; respectively, have full rank p. Due to 
the invertibility assumption on 0, it follows that these domains are residual subsets (countable 
intersections of open dense subsets) of C pxp (z). Note that a linear fractional transformation is 
fully determined if it is specified on such a residual set. 

It is easily established that the following group properties hold for the LFTs associated with 
two invertible matrices $ and 

7i> o 7* = 7i»i», (39) 
7*o 7* = 7*$. (40) 

Here it may be noted that the compositions of mappings on the left hand sides of these equations 
may formally have domains which are strict subsets of the domains of the mappings on the right 
hand sides. However, since the domains of the compositions of mappings on the left-hand sides 
are residual sets too, this issue is resolved by identifying (compositions of) LFTs if they agree 
on the intersection of their domains, which will be our policy throughout this paper. 

It then holds that the LFTs 7e and 7e are bijections from their domains A4@ and .Me to their 
co-domains .M@-i and A^e-i, and that every LFT can be represented in each of the two forms 
described above. We have the following lemma. 

Lemma 4.1 Let G C 2px2p (z) be an invertible rational matrix. Then: 

(i) 7e and To are bijections. Their inverses are given by 

Tq 1 = Tq-i = Tjej, 
Tq 1 = Tq-i = Tjbj- 

(ii) 7e = fj e -ij and fe = T je -ij. 

We now proceed to study the conditions under which two LFTs coincide. From the literature, 
the following result is well known |26j . 

Lemma 4.2 Let <E>,\I/ G C 2px2p (z) be two invertible rational matrices. It holds that 7$ = 7* if 
and only if there exists a scalar function A G C(z) for which $ = 

For our purposes, however, it will be convenient to develop a slightly specialized version of this 
lemma, which states that the same result remains to hold true if two LFTs 7$ and Tq/ are merely 
required to coincide on the subset of p x p lossless functions. 

Lemma 4.3 Let <I>,\E' £ C 2px2p (z) be two invertible rational matrices. If the p x p lossless 
functions are all contained in M.$ f] M.^, and if for all lossless G it holds that 7$(G) = T^{G), 
then there exists a scalar function A G C(z) for which <3? = A^. 
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Proof. Note that upon application of the bijection 7yp-i, we have that the identity 7$(G) = 7*(G) 
becomes equivalent to 7^,-i$(G) = G. Denoting = ^L'~ 1 <I>, this can be rewritten equivalently 
as 4 G + 3 = G0 2 G + G@i, since (0 2 G + ©i) -1 exists (for all lossless G). If G is lossless, 
then also — G is lossless. Comparing the two expressions which follow from substitution of these 
values for G, it is found that for all lossless G 

9 4 G = GGi, 
G 3 = G0 2 G. 

Substitution of G = I p yields @i = 04 and ©2 = ©3. Setting G = which is again lossless, 

yields 3 = z~ 2 Q 2 - In combination with ©2 = ©3 it follows that ©2 = ©3 = 0. 

We are left with the identity @iG = G@i for all lossless G. Note that every unitary matrix is 
lossless, so that G can be varied over all the products of signature matrices and permutation 
matrices. From this it easily follows that ©1 is a scalar multiple of the identity matrix. So there 
exists a scalar function A such that Q% = ©4 = XI p . (Rationality of A is obvious from rationality 
of ©, which is a consequence of rationality of and ^.) Hence = Xl2 P , and therefore $ = A^. 
□ 

Remark. The proof of Lemma f4.3l is such that it also applies to a setting in which all the rational 
matrix functions are restricted to be real. In the complex case, one could replace the choice 
G = z~ 1 I p (which is lossless of McMillan degree p) by the choice G = il p , so that the whole 
proof involves values for G from the set of unitary matrices only (i.e., G lossless of McMillan 
degree 0). In the real case with p > 1 a more careful analysis shows that a proof can also be 
designed using orthogonal matrices only. However, in the real case with p = 1, the only lossless 
functions of McMillan degree are G = ±1. Here, lossless functions of McMillan degree > 1 are 
needed to obtain the desired result. 

For our purposes, it is important to study in some more detail the situation where the 2p x 2p 
rational matrix is J-inner and the p x p rational matrices G on which the LFT 7e acts are 
lossless. 

Proposition 4.4 Let be J-inner of size 2p x 2p and of McMillan degree m. If G is p x p 
lossless and of McMillan degree n, then G S A^e an d the matrix function G = 7e(G) is also 
lossless and of McMillan degree <n + m. 

Proof. The proof can be easily adapted from [31 sect. 18. 2.]. A realization of G with (n + Tri- 
dimensional state-space can be easily computed (see formula (4.85) in [IS]), which proves that 
the McMillan degree of G is at most n + m. □ 

Also, if the LFT corresponds to a J-inner matrix, we have the following elementary results. 

Lemma 4.5 Let Q(z) be J-inner of size 2p x 2p. 
(a) The mappings 7e andT®t coincide. 
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(b) If P and Q are p x p unitary matrices, then @p t Q(z) := @(z) 
and for all G £ Mq pq : 

T QpMz) {G{z)) =T e(z) (QG(z)P*). 



P 
Q 



is again J -inner, 



(c) Let K 



I p 
L h 



Then the function defined by 



&°{z) := KQ(z~ 1 )K 
is again J -inner, and for all G 6 Mq it holds that T@(G)* = T@°(G* 



(41) 



Proof, (a) Since is J-inner, it holds that 7ett = Tj®-ij = Te- 

(b) The proof that Qp^q(z) is again J-inner is elementary. The rest of this statement follows 
from the group property of LFTs and from the simple action of an LFT associated with a block- 
diagonal matrix. 

(c) Note that Te(G)* = (G*Q* + e;)" 1 ^^ + G*) = Tkq*k{G*) = Tj Km -i KJ {G>) = 
Tk(q*)»k(G*) = 7e°(G*). The proof that 0°(z) is again J-inner is elementary upon observing 
that KJ = -JK. □ 



A particular situation of interest occurs when the J-inner matrix function Q(z) happens to be of 
McMillan degree 0. In this case &(z) = M is a constant J-unitary matrix and the associated class 
of LFTs is the one of generalized Mobius transformations Tm- The following theorem indicates 
how the action of such transformations Tm on lossless functions G(z) can be represented entirely 
in terms of balanced state-space realizations. 



Theorem 4.6 The linear fractional transformation Tm associated with a constant J-unitary 
matrix M is a bisection on the set of lossless functions which preserves the McMillan degree. 
Mi M 2 
M 3 M 4 

minimal state-space realization (A, B,C,D). Then, a minimal state-space realization (A, B,C, D) 
for G(z) = Tm(G(z)) is given by: 



Let M 



be the block decomposition of M and let G(z) be a lossless function with 



A 
B 

C 

b 



A - B(M 2 D + Mi)- 1 M 2 C, 
B(M 2 D + M 1 )-\ 

[M 4 - (M 4 D + M 3 )(M 2 D + M 1 )- 1 M 2 ]C, 



(42) 



(M 4 D + M 3 )(M 2 D + Mi)" 1 . 
If in addition (A,B,C,D) is balanced, then this realization (A,B,C,D) is also balanced. 



Proof. A J-unitary constant matrix M is in particular J-inner and its inverse is also J-inner, so 
that Tm is clearly, by Prop. 14.41 a bijection on the set of lossless functions. 
The computation on the realization is just a particular case of formula (4.85) in [18]. This shows 
that the McMillan degree does not increase under the action of Tm- It cannot decrease either, 
because the inverse mapping Tm- 1 = T M l can be treated likewise, and also does not increase 
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the McMillan degree. Hence, the McMillan degree is preserved by 7m- 

Finally, recall that a minimal realization of a lossless function is balanced if and only if the asso- 
ciated realization matrix is unitary. Starting from a minimal balanced realization (A, B, C, D) it 
now can be verified by direct computation that the realization matrix associated with (A, B, C, D) 
is also unitary and thus the realization is balanced. □ 



5 The tangential Schur algorithm 



In this section we outline the use of the tangential Schur algorithm for the recursive construction 
of a parametrization of the space of stable all-pass systems of fixed finite order. It is derived from 
the method of [UHO], where the tangential Schur algorithm is used to construct an infinite atlas 
of generic overlapping parametrizations for the space of p x p inner functions of McMillan degree 
n. The relationship between these two situations is constituted by the map R(z) i— > R(z)^ 1 , 
or equivalently R{z) \-t R?(z) (see also Eqns. (fTTj) and (fT7|) ). which is used to relate the space 
of p x p inner functions of McMillan degree n to the space of p x p stable all-pass systems of 
McMillan degree n. 



In the context of inner functions, the tangential Schur algorithm consists of an iterative procedure 
by which a given p x p inner function of McMillan degree n is reduced in n iteration steps to a 
p x p inner function of McMillan degree (i.e., to a constant unitary matrix). In each iteration 
step the McMillan degree of the inner function at hand is reduced by 1, by application of 
a suitable linear fractional transformation which is chosen to meet a particular interpolation 
condition, and which involves an associated J-inner matrix function of McMillan degree 1 (or 
elementary J-inner factor, see Section [2]). The actual parametrization procedure consists of 
the reverse process, by which a chart of inner functions of McMillan degree n is constructed 
in n iteration steps, starting from an initial unitary matrix. The choice of interpolation points 
and (normalized) direction vectors may serve to index such a chart, while the local coordinates 
correspond to the n Schur vectors. 



In the Schur algorithm, the elementary J-inner factors involved are of the form (|28|) and it 
suffices to consider those which are analytic on the closed unit disk (see the remark at the end 
of the section). Such an elementary J-inner factor Q(z), having its pole outside the closed unit 
disk at z = l/w, can be represented as: 



@(u,v,w,£,H)(z) 



h P + 



b w {z) 

MO 



u 




u 


* 








J 


V 




V 




(1 




IMI 2 


) 



A 



(43) 



where 

(1) w G C, \w\ < 1, 

(2) u e C pxl , (Ml = l, 

(3) v e C7 xl , |H| < l, 
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(4) £ e C, |e| = l, 

(5) H is a 2p x 2p constant J-unitary matrix. 



This representation is obtained from (|28p in which some normalizing conditions have been im- 
posed: 



(i) the vector x has been rewritten x = a 
thus Hull < 1, since (|29j) must be satisfied 



u 
v 



, where a is a nonzero scalar, ||u|| = 1 and 



(ii) an arbitrary complex number £ on the unit circle has been introduced which allows to 
write the right J-unitary factor H as H = 0(£). 



Since any nonzero scaling of the vector x does not affect the expression (128|) . we have that 

@(u, v, w, £, H)(z) = X x (Mz)MCr 1 ) H = X X (b w (z)) X x (MO)" 1 H, (44) 

and @(u,v,w,£,H)(z) is of the form (j28|) up to a right constant J-unitary factor. The constant 
matrix X x (b w (£) ) = Xe(&uj(£)) _1 i s indeed J-unitary, since |£| = 1 yields 1 6^(^)1 = 1 ( see 
Lemma 12.211 . 



Once £ is chosen arbitrarily on the unit circle, the representation (143p is unique up to a uni- 
modular complex number, which corresponds to the remaining freedom in the scaling of the 
vector x. It will turn out that for several of our purposes, the additional freedom to choose £ 
and H is of rather limited use. In the literature on the parametrization of inner functions it 
therefore often happens that fixed choices are made, in particular £ = 1 and H = I2 P (see, e.g., 
[TJ[T0]). However, especially the role played by H cannot be ignored in establishing the precise 
connection between the tangential Schur algorithm and the state-space approach involving the 
mappings J-jjy i n Section [6j Also, for the construction of atlases of overlapping parametriza- 
tions, the additional freedom to choose H leads to useful new alternatives, as will be discussed 
in Section [3 Therefore, we will not fix a choice yet and present our results for general £ and H. 

The following lemma establishes a useful new factorization result for the class of J-unitary 
matrices X x {a) with x*Jx > 0. 



Lemma 5.1 Let x G C 2pxl be such that x*Jx > 0. Let x be written as x = a 



u 

V 



, where a is 



any nonzero scalar such that u, v € C pxl satisfy \\u\\ = 1 and \\v\\ < 1. Lt then holds that 



X x (a) = H(uv*)X 



u 




(a)H(uv*y\ 



(45) 



where H(uv*) denotes the Halraos extension of the strictly contractive matrix uv* . 



Proof. Note that a, u and v with the required properties always exist, making uv* into a strictly 
contractive matrix. The Halmos extension H(uv*) is not difficult to compute explicitly as the 
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Hermitian J-unitary matrix 



H(uv*) 



)UW 



=uv 



v*=ii 



=t)M 



Ip - (1 ! 



v*=ii 



\ vv 
• IMI 



(46) 



and it satisfies 



-v 



H(uv*) 

It then follows that X x (a)H{uv*) = H(uv*) + 



a/i - IMP 

(a-l) 







(«-l) 
proves the lemma. 



It 














H(uv*) + (a-l)H(uv 



(i - Ik! 

it 











1i 









u 




JH(uv*) = H(uv*) + 
J /ICuu*)X r t (a), which 







□ 



The following proposition plays a central role in our construction of balanced parametrizations 
for discrete-time lossless transfer functions. 



Proposition 5.2 The J -inner matrix function @(u,v,w,£, H)(z) can be factorized as: 

Q(u,v, w,£,H)(z) = H(uv*)S u , w (z)S u , w (b w {t))~ l H{uv*)- l H, (47) 
where H(uv*) denotes the Halmos extension of the strictly contractive matrix uv* , and where 



Su,w(z) '• — X 



u 




I p - (1 - b w (z)) uu* 
/„ 



(48) 



Proof. This is an immediate consequence of the expression (|44p of &{u, v, w, £, H)(z), in combi- 
nation with Lemma l5.lt in which a is replaced by the Blaschke factor b w (z). □ 

The next proposition establishes that application of the linear fractional transformation associ- 
ated with the elementary J-inner function 0(u, v, w, £, H)(z) to a lossless function of McMillan 
degree n, results in a lossless function of McMillan degree n + 1 which satisfies an interpolation 
condition involving only u, v and w. 



Proposition 5.3 Let Q(z) = 0(u, v, w, £, H){z) be the elementary J -inner function defined in 
Eqn. |^3| ), with u, v, w, £ and H satisfying the accompanying properties. Let G{z) be a lossless 
function of McMillan degree n. Then 



G(z)=T e{z) (G(z))=T @Kz) (G(z)) 
is lossless of McMillan degree n + l and satisfies the interpolation condition 

G(l/w)u = v. 



(49) 
(50) 
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Proof. Proposition 14.41 asserts that G(z) is lossless. From Proposition 15.21 we have that 



©(*) 



° 7m, 



where M = Su^Q -1 H (uv*)~ l H . The linear fractional transformations associated with the 
J-unitary matrices M and H(uv*) are both generalized Mobius transformations which leave 
the McMillan degree unchanged, see Theorem 14.61 Only the linear fractional transformation 
associated with the matrix S UjW (z) can change the McMillan degree, but it has a simple form: 
for the lossless function R(z) := Tm(G(z)) of McMillan degree n it holds that 

■= r SuMz) (R(z)) = R(z)(i P - (i - bl(z)) uu *). 

Because I p — (1 — bw{z))uu* is the transfer matrix of a lossless system and, by ([2T]) . det R(z) = 
b w {z) det R(z), it follows from Lemma 12. II and from the fact that stability ensures that no com- 
mon factors can occur, that R(z) has McMillan degree n + 1. Then also G(z) := Tjjf uv *\(R(z)) 
has McMillan degree n + 1. 



For the matrix X x (a) it holds that x*JX x (a) = ax*, J. Setting x 
this results in the following identity, in view of Eqn. (|44|) : 

[ u* -v* } e(w) = o. 

Application of the operator " yields the equivalent relationship 



u 
v 



, a = b w (z) and z 



w 



u 
-v 



0. 



(51) 



(52) 



Note that G(z) = Tq%i z \(G[z)) satisfies the relation 

G{z) I p ] = (G(z)Ql(z) + Qi(z))- 1 [ G(z) I p ] e\z). 



(53) 



The matrix G(z)@ 3 (z) + ®\{z) is analytic and invertible everywhere outside the open unit disk, 
which follows in a manner entirely analogous to the proof of invertibility of @2(z)G(z) + Q\(z) 
at points of analyticity outside the unit disk, see the proof of Proposition 14.41 Setting z = 1/w 
the interpolation condition f)50f) now follows. □ 



Remark. Using relation (|5ip and the fact that G(z) = T@( z ){G{z)) can be rewritten in the form 



G 



e 



G 



(e 2 G + e 1 )-\ 



yield an alternative form of the interpolation condition (1501) : 

u* = v *G(w) 



(54) 



(55) 



provided that the matrix Q2(w)G(w) + @i(w) is invertible. Equivalently, this requires that w 
does not show up as a pole of G. Only then the form (|55l) is valid, while the form ([50]) applies 
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always. This illustrates the usefulness of the identity 7@tt = 7e when dealing with analyticity. 
In a similar vein, note that for w = the function @(u,v,w,£, H) is not proper and does not 
admit a state-space realization. However, the matrix function &(u, v, w, £, does admit a 
state-space realization, and formula (4.85) in [18] can be used to compute the action of the 
associated LFT entirely in terms of realizations. 

The following proposition addresses the reverse process and constructs from a lossless matrix 
function G(z) of McMillan degree n + 1, a new lossless matrix function of reduced degree n. To 
achieve this, it proceeds from an interpolation condition of the form (|50p . Though this result is 
well-known (see, e.g. [3], pQ), we find it of interest to give a more constructive proof of it. 



Proposition 5.4 Let G be a px p lossless function of McMillan degree n + 1 which satisfies an 
interpolation condition of the form 

G(l/w)u = v, (56) 

in which w G C is an interpolation point with \w\ < 1, u G C pxl is a direction vector with 
\\u\\ = 1 and v G C pxl is a Schur vector satisfying \\v\\ < 1. 

Let £ G C be an arbitrary number of modulus 1, and let H be an arbitrary constant J -unitary 
matrix. Then G admits the representation 



G — Te( u ,v,w,£,H){G)) 
for some p x p lossless function G of McMillan degree n. 



(57) 



Proof. The interpolation condition (|56|) can be rewritten as follows: 
G{l/w)u = v O 0=\d(l/w) L p 



u 
-v 



^ 







G(l/w) I p 



H(uv*)- l H(uv*) 



G(l/w) L p H(uv*) 



u 




& = T H{uv * ) -i(G)(l/w)u. 



u 
-v 



The lossless function R := Th( U v*)-' l {G) = Th( U v*)- 1 (G) has the same McMillan degree n + 1 as 
G and satisfies the interpolation condition 

R(l/w)u = 0. 

This means that R(z) has a zero at z = 1/w, or equivalently a pole at z = w. Now let B w (z) 
be defined as 

B w {z) := L p - (1 - bi(z))uu* = (L p - (1 - b w (z))uu*r\ (58) 

which is p x p lossless of McMillan degree 1, satisfying det B w (z) = bw(z) and having its pole at 
z = w. Then R(z) can be factored as 



R(z) = R(z)B w (z), 



(59) 
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where by construction R{z) is p x p lossless of McMillan degree n. (This is a matrix version of 
the Schwartz lemma on which the Potapov factorization of lossless and inner functions is based; 
see [23].) The identity (f59l) can be rewritten as 

R = T Su!W (R) 

so that finally 

G = Te( u ,v,w,z,H){G), 

where G = Tm-i(R), with M = S^^) -1 ^ (-u-u*) _1 -ff, is also lossless and of McMillan degree 
n. □ 

Remark. In the standard case with w = the value of G(oo) corresponds to the direct feedthrough 
term D of any state-space realization (A,B,C, D) of G(z), so that the interpolation condition 
then takes the form Du = v. 

The tangential Schur algorithm can now be described. It is used in Section[7]to construct atlases 
of generic charts for the manifold of p x p lossless systems of degree n. 



Theorem 5.5 (Tangential Schur algorithm) Let G^ be a p x p lossless transfer matrix of 
McMillan degree n. For k = n, . . . , 1, let interpolation points w k £ C be given with \w k \ < 1, 
let constants G C be given with |^| = 1, and let mappings : (u,v,w,£) i-> Hf.(u,v,w,£) 
be given which assign a constant 2p x 2p J -unitary matrix Hk(u,v,w,£,) to each quadruple 
(u,v,w,ti) ePxPxCxC with \\u\\ = 1, \\v\\ < 1, |w| < 1 and \£\ = 1. 
Then for k = n, . . . , 1, there exist vectors Uk £ C p with \\uk\\ = 1, such that the vectors Vk £ C p 
constructed recursively by the following formulas, all have length \\vk\\ < 1." 

v k :=G^(l/w k )u k , (60) 
Q k := @{u k , v k , w k ,i k , H k (u k , v k , w k , £ k )), (61) 
Q(k-l) :=Tq^(G^). (62) 

With such a choice of the unit vectors u k (k = n, . . . , 1) each of the functions G^ is lossless of 
McMillan degree k and one can write 

<?W = Te n (Te n .A- ■ ■ (TeAG {0) )) ...)) = Te^-eAG^), (63) 
where G^ is a constant unitary matrix. 



Proof. In view of Proposition 15. 4\ which is applied repeatedly, it only remains to show that if G 
is a lossless function of degree n > 0, then for every choice of interpolation point w with \ w\ < 1, 
a direction vector u with ||n|| = 1 exists for which the Schur vector v := G(l/w)u satisfies 
\\v\\ < 1. A function theoretic proof of this fact is given in [X0^ Lemma 6]. An alternative 
state-space argument runs as follows. 

For < 1, the mapping z i — y b— w (z^j constitutes an automorphism of the unit disk. Thus, 
the function G(z) := G(b- W (z)) is again lossless (and also of McMillan degree n > 0). This 
shows that G(l/w) = G(oo) appears as the direct feedthrough matrix D of any state-space 
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realization (A, B, C, D) of the lossless function G. If such a realization is chosen to be minimal 
and balanced, then D is a proper submatrix of the associated unitary realization matrix; see 
Proposition 13.21 Clearly, some column of D must have length < 1, since otherwise all columns 
of B would be zero, contradicting minimality with order n > 0. Thus, a suitable vector u can 
always be found among the set of standard basis vectors {e%, . . . , e p }. As a matter of fact, any 
normalized basis will do. □ 

Remark. Similar results can be proved in which the interpolation condition (150j) . which involves 
a relationship among the columns of G, is replaced by an interpolation condition involving its 
rows. This can be done by applying Proposition 15.41 to the lossless function G*(z) and using 
part (c) of Lemma 14.51 The interpolation condition then takes the form 

u*G{l/w) = v*, (64) 

and should not be confused with the alternative form (|55p for the interpolation condition (|50p . 
As a counterpart to the result of Proposition [5741 one also obtains an alternative linear fractional 
representation 

G = Tq( u ,v,w^H)°{G) 

for some pxp lossless function G of degree n. Observe that the J-inner function Q(u, v, w, £, H)° 
is now analytic outside the open unit disk. Thus, the theory of this section admits a dual theory 
in which elementary J-inner factors analytic outside the disk and interpolation conditions on 
the rows are used. 

Note that the tangential Schur algorithm described above proceeds from a set of n interpolation 
conditions on the columns of a sequence of intermediate lossless functions. Obviously, a set 
of n interpolation conditions on the rows of a sequence of intermediate lossless functions can 
also be imposed. But since each iteration step involves a single interpolation condition only, 
mixed situations can also be treated in which n interpolation conditions are prescribed involving 
columns and rows. 



6 Connection between the classes of mappings Tq( z \ and T\jy 

In the scalar case, mappings of the form T\jy with U = I2 have been used in [T5] to recursively 
construct a balanced canonical form for the space of discrete-time lossless systems of finite 
McMillan degree. The parameters that occur in this recursion have the interpretation of Schur 
parameters, corresponding to the situation with interpolation points w at zero. With this 
connection in mind, it is the purpose of the present section to clarify the relationship between 
the two classes of mappings Tq( z ) and J-\jy, with O(z) a J-inner matrix function of McMillan 
degree 1 and U and V unitary. We are interested in investigating the possibilities for representing 
a mapping 7e( 2 ) m terms of a corresponding mapping J-jjy- This will give us balanced state- 
space parametrizations directly in terms of the set of parameters u, v, w, £ and H used in the 
tangential Schur algorithm. Moreover, unitary matrices will be involved in the computation of 
these realizations, and these are known to be numerically well-conditioned. 
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Let Cn denote the set of rational p x p lossless transfer matrices of McMillan degree n and let 
= 



Theorem 6.1 Let U and V be unitary (p + 1) x (p + 1) matrices, block-partitioned as 



U 



ku ft 



V 



ay My 

ky f3*y 



(65) 



with kjj and ky scalar and the other blocks of compatible sizes. 

The mapping Tu,v : C p — > C p can be expressed as a linear fractional transformation 7$ : C p — > 
C p for some J -inner function $ if and only if (ku, ky) ^ (0, 0). 

(i) If\kjj\ = \ky\ = 1, then Mjj and My are unitary matrices and J-jjv ^ s generalized Mobius 
transform Tuy(G) = MyGM^ which can be represented as 7~$(G) with a constant J -unitary 
matrix 

(ii) If (kjj,ky) 7^ (0,0) and \kjj\ < 1 or \ky\ < 1, then J-jjv can be represented as 7$ with an 
elementary J -inner function <3? having a pole at ky /ku, located at infinity for ku = 0. 

Proof. Because U and V are unitary, according to Proposition 13,41 the mapping Tuy maps C p 
into C p . If $ is J-inner, then the mapping 7# also maps C p into C p but in addition it is bijective, 
see Lemma [4. II and Proposition 14.41 Thus, a necessary requirement for Tuy to be representable 
as a linear fractional transformation 7$ with a J-inner function <1>, is injectivity. 
First consider the scalar case p = 1 with (ku,ky) = (0,0). Then U and V are 2x2 diagonal 
matrices and all lossless G are mapped to the unimodular constant ayo^j. Therefore the mapping 
Tuy obviously is not injective. 

In the multivariable case p > 1 with (ku,ky) = (0,0), it also holds that Tuy is not injective 
and thus cannot be represented as an LFT associated with a J-inner function. Indeed, let 
P(0 = Ip + (£ - tyPuPu, and = I p + (£ — l)/3 v f3 v where £ / 1 is a unimodular number. 
Then and Q(£) are unitary and we have 



U 



so that by Lemma 



1 
P(0 



I P 



U, V 



1 

Q(0 



o £ 

(i), (iii), for every lossless function G, 



I P 

o e 



Tuy{Q{i)G{z)P{iY) = T 



' Ip 


" 






" 










£ _ 



(G(z)) = Tuy(G(z)). 



Therefore injectivity of Tuy requires that Q(£)GP(£)* = G for every lossless function G. How- 
ever if this is the case, then (3yf3 v G = Gfiufiu ^ or an lossless G. Choosing G = I p this implies 
PvPy = fiufiu- Choosing G = I p — 277*, where 7 G C p is a vector of length one which is 
independent of f3y and not orthogonal to j3y, it follows that /3y/3y77* = 77*/3[//3^}, which yields 
a contradiction. Therefore Jj/y is not injective. 
Next, assume that (ku,ky) 7^ (0,0). Let 



Ji 



Ip 




and J2 




1 
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and consider the linear fractional transformation C defined by 

C(H(z)) = (Jitf (*) + J 2 )(J2H(z) + Ji)- 1 . 



This linear fractional transformation is well-defined for every H(z) such that, in the block 
decomposition of H(z) 

' H^z) H 2 {z) 
H 3 (z) H A {z) 

with Hi(z) scalar, H^z) is not identically zero. It can be computed as 



H(z) 



C(H) 



H\ — H 2 H A 1 H% H 2 H A 



H, 



-i 



'4 ^6 "4 

The mapping Tyy can now be written using this transformation as: 



F uy {G{z)) = [l p Q]jr[V 



1 

G{z) 



U* -z 




1 



Now, 
£ IV 



1 

G{z) 

1 

G(z 

ay MyG 
ku P*u 



U* -z 



Ji(V 




1 


1 



U )+ J 2 U 



h v 



1 

G(z) 




1 



U)+J 1 U 



a v Mu 
ky - zku /3 V G - z/3 v 



Using the following expression for the inverse of a partitioned matrix 



Mi M 2 
M 3 M 4 



M 3 




-1 1 



+ 



-M7 L M A 



(M 2 - MiM 3 " 1 M 4 )' 1 [ I p -MiM 



-i 



in which M 2 and M3 are square submatrices, and M 3 1 is assumed to exist (which holds by 
assumption in our case), we obtain 



Fuy{G)=[ay MyG 

ay (pyG-z^ 



fj v G-zfj* v 

ky—Z ku 



M v - 



au (PyG-zPb) 



My G - 
My - 



ky — zku 
ay P v 



G + 



M v - 



za v p*j 



ky — zku J ky — zku 

And finally, Fu,v{G) = 7$(G) for $ given by 

zaufiu 



au (PyG-zPtj) 
ky — zku 

au Py 
ky — zku 



ky — zku 
-1 



G + ( Mu + 



zau Pjj 
ky — zku 



-1 



M v + 



ky — zku 
zayf3 v 

ky — zku 



My - 



ky — zku 
ayfiy 



ky — zku 



(66) 
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Exploiting unitarity of U and V, it then is straightforward to establish that 

J - J$(A)* = __ 

(fcy — 2 ku)(kv — A fej/) 

from which it follows that &(z) is J-inner. The expression (|66|) for 3>(z) is of McMillan degree 
1 having a pole at ky/kjj (located at infinity if kjj = 0), except if \kjj\ = \ky | = 1, in which case 

M v 



au 




au 


ay 




ay 



$(z) is of McMillan degree since au = ay = and <&(z) 
□ 



My 



becomes constant. 



The previous theorem makes clear under which conditions a mapping J~uy with unitary matrices 
U and V admits a representation as an LFT associated with an elementary J-inner factor 
Depending on the location of the pole ky/ku, this factor takes the different forms described in 
Theorem 12.31 In the cases where such a representation exists the proof is constructive, yielding 
the possible expression (166 p for <£. Using Lemma [4.3l all the possibilities for $ can be computed. 



Lemma 6.2 Let 3>(z) and ^(z) be elementary J-inner factors for which 7$ and 7\p agree on 
C p . If p > 1 then there exists a unimodular constant p such that 

<&(z) = p^{z). 

If p = 1, then either there exists a unimodular constant p such that &(z) = p^{z), or &(z) has 
a pole at z = w with \w\ ^ 1, *&(z) has a pole at z = 1/w and there exists a unimodular constant 
p such that 

= pbt(z)*(z). 

Proof. By Lemma l4.3| 

= b{z)V{z), (67) 

for some scalar function b(z), and since <3? and ^ are J-inner, b{z) must be a Blaschke factor. 
Computing the determinants gives 

det*(z) = b{z) 2p det^(z). 

But the determinant of an elementary J-inner function is either equal to a unimodular number 
(if the pole is located on the unit circle) or to a Blaschke factor, having McMillan degree 1. If 
p > 1 it therefore follows that that b(z) must be a unimodular constant. If p = 1 an alternative 
possibility occurs if $ has a pole at z = w with \w\ ^ 1. In that case ^ may have a pole at 
z = 1/w and b(z) may also be chosen equal to pb\ iJ (z) for some complex number p of modulus 
one. If ^ is of the form 

= (I 2 - (1 - b w (z))x(x*Jx)~ 1 x*J) H, 

then 

= p (/ 2 - (1 - U(z)My* Jy)- 1 ^ J) F, 
where y is a solution to y* Jx = 0. □ 

As a consequence of Theorem 16.11 and Lemma 16.21 we have the following result. 
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Corollary 6.3 Let U and V be unitary matrices of size x (p+1), partitioned as in 1165]) . 

If p > 1, then the mapping T\jy can he represented in the form of a mapping Tq( u ,v,w^,h) if and 
only if\kjj\ < \ky\. Likewise it can be represented in the form of a mapping Tq{ u ,v,w^,h)° if an d 
only if | ku | > | ky \ ■ 

If p = 1, a representation of J-jjy *n the form of a mapping Tq( u ,v,w,(,,h) exists if and only 
if a representation in the form of a mapping 7q( u vw £, H)° exists, which holds if and only if 



Proof. The statements for the case p > 1 follow from the earlier results as indicated. For p = 1, 
note that T\jy = J~vjj which follows from Lemma [3.3l (ii) upon complex conjugation and noting 
that transposition has no effect on a matrix of size lxl. □ 

To complete the picture it remains to characterize the elementary J-inner factor analytic 
inside the disk in case \kjj\ < \ky\. As we have seen, if a mapping J-jjy is given, the expression 
(|66p together with Lemma 16.21 provide all the possible alternatives for <£. Conversely, if an 
elementary J-inner factor $ analytic inside the disk is given, one may wonder whether 7$ 
admits a representation as a mapping of the form T\jy. Here it will prove to be essential to 
introduce the J-inner matrix function Q(u,v,w)(z) defined by 

e(u,v,w)(z) = H(uv*)S u>w {z)H(wuv*). (68) 

The main result of this section can now be stated as follows. 



Theorem 6.4 Let Q(z) be an elementary J-inner factor analytic inside the unit disk. 7/7e( 2 ) 
coincides with a mapping Tyy with U and V unitary, then Q(z) can be written as 



00) = Q(u,v,w)(z) 



P 
Q 



(69) 



for some u, v G C p with \\u\\ = 1, \\v\\ < 1, some w € C with \w\ < 1, and some p x p unitary 
matrices P and Q. 



In that case one can take U = U 



1 
P 



and V = V 



1 
Q 



, where 



U 



V 



v^-H 2 u j fi i W VHMI 2 
^i-H 2 |M| 2 ' p ^/i-h 2 |M| : 



^i-H 2 IM! 2 

Jl- |ui| 2 ||ti||' 



— \w\ 2 \\v\\ 2 



it. 



v I P -(1 



\/i-h\ 2 \M 2 R 



(70) 



(71) 



Proof. Theorem 16.11 ensures that a mapping J-jjy can be represented by a linear fractional 
transformation for the J-inner elementary factor 3>(z) given by (j66f) . If p > 1, it follows from 
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Lemma E21 that \ku\ < \ky\- If P = 1, since the mappings Tjjy and Tyjj coincide, we still may 
assume \kjj\ < \ky\, and in both cases we have that 

$(z) = pQ(z), 

for some unimodular number p. 

Since \kjj\ < \ky\, fcy and k\j ^ 1, so that \\ajj\\ 7^ 0. Let 



w - — 

ky ' 



ay 
I «t/ 1 



ay 



(72) 



which satisfy \w\ < 1, ||u|| = 1 and \\v\\ < 1. The function Q(z) is, by Proposition 15. 2\ of the 
form 

Q(z) = H(uv*)S UjW (z)H(E) £ ° 

for some strictly contractive matrix E and some unitary matrices P and Q, depending on the 
matrices U and V. To see this, note that if Q(z) is an elementary J-inner factor written in the 
form @(u,v,w,£, H)(z) having its only zero at z = w, then the matrix Q(w) has a nontrivial 
left kernel of dimension one, which is characterized by Eqn. (15ip . For the matrix function $>(z) 
given by (j66[) it is easily shown that [ u* —v* ] 3>(w) = for w, u and v defined by Eqn. (|72p . 
Note that, the matrices U and V being unitary, 

\\au\\ 2 + \w\ 2 \ky\ 2 = 1 
l|at/|| 2 |bl| 2 + \ k v\ 2 = 1, 



so that 



at/ 



u/ 



|Av| 



0~ 



\w\ 2 \\v\\ 2 



(73) 



yl — |u>| 2 ||t;|| 2 a/1 

Now, the quickest way to proceed is to compare the realizations of 3?" and 0^ (note that will 
fail to have a realization if w = 0). Using (j66[) we get 



M* 







Mi 



+ 



C *A-1 



Am 



while a state-space realization of Sl :W (z) is easily computed to yield that 

S\z) =V + C(z-w)- 1 B, 



ky 



ky 



where 



C 



V 

P* 
Q* 



P* 
Q* 



H(E) 



In 



(1 + w)uu* 








H(uv*), 



H(E) 



\w\ 2 u 







and B=[ y/l - \w\ 2 u* ] H(uv* 



Comparing the (l,2)-block-entries of the expressions for the direct feedthrough yields, using 

(dD, 

P*(I P - EE*y 1/2 (E - wuv*)(I p - w*y 1/2 = 
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and thus E = wuv* or equivalently 

@(z) = e{u,v,w)(z) 
as claimed. Replacing E by its value gives 



P 
Q 



V 



P* 
Q* 



1 + 



i-M 2 



uu 



C 



to 



y/l — \w\ 2 \\v || 2 





MlV 1 




U! 2 ||?j|| 2 









It 





Q* 




WV 


riven 


by 







vu 



Ip- 1 1 





^/i-M 2 IMI 2 



and B = ^ S [ u* 



The matrix Mjj is therefore given by 



M v = p \I P - 11 + 



\J\ — |u)| 2 ||t>| 12 



uu* P. 



From (1721) and (1751) we have that 



\kv\_ 

ky 



which determines the transformation between the two realizations, so that 

ky 



U k v Pv 



P 



\k V \ 



C* = p 



0" 



\w\ 



\k v \ Vl-\w\ 2 \ 



[ u* 



W V 



Finally it is obtained that the matrices U and V are of the form 



U 



I P _o 
fe 



U 



1 

pP 



and V 



I P 
u M 



P 
Q 



V 



1 

pQ 



with [/ and V specified by (|70|) and ([7T]). Since a pre-multiplication by 



I P 
fe 



1 

pip 

the values of U and V given in the theorem are obtained 



post-multiplication by 



and a 



do not change the mapping T\j y according to Lemma | 



□ 



7 Overlapping canonical forms for lossless and stable systems 

In the tangential Schur algorithm as presented in Theorem 1 5. 5 ( a lossless function G = of 
McMillan degree n is broken down to a constant unitary matrix G^ using a sequence of inter- 
polation points Wk, direction vectors Uk, Schur vectors Vk, unimodular constants and constant 
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J— unitary matrices Ufc, Wk,Ck), for k = 1, . . . , n. The set of values for w^, Uk, S,k and the 

set of mappings at the n recursion steps can serve to index a generic chart for the differen- 
tiable manifold of lossless functions of degree n, if the mappings are sufficiently smooth. The 
Schur parameter vectors v& together with the unitary matrix that we finally reach, then 
provide the local coordinates for this chart. An infinite atlas of overlapping generic charts is 
obtained by varying the choices for the u^, w^, and H^. Sub-atlases can be extracted from it, 
as the one built in pQ which employs elementary J-inner matrices of the form 0(iik,Vk,Wk, l,I p )- 
Such a choice implies that G^ is the value of the lossless function G at 1. We may even take 
a finite atlas by fixing all the interpolation points at zero and letting the values of the Uk vary 
among a canonical basis of C p (cf. the proof of Theorem \5.5\i . However it is interesting to keep 
in mind the richness of the possibilities, including the one to mix interpolation conditions on the 
columns (I50p with interpolation conditions on the rows (|64p . This could be of interest in order 
to describe lossless transfer functions with a particular structure. 

Now let constants be chosen with = 1 for all k = l,...,n, and let the mappings 

H k (u,v,w,£) = H(uv*)S U)W (S,)H(wuv*) 

be chosen so that &(u,v,w,^k,Hk(u,v,w,^k)) = @(u,v,w) for all u G C p , ||u|| = 1, v G C p , 
\\v\\ < 1, w G C, \w\ < 1. Then one obtains an atlas for the manifold of lossless systems that 
can be described directly in terms of state-space realizations using the T\jy mappings. 

Let us now describe such an atlas for the manifold Cn of p x p lossless functions of McMillan 
degree n in more detail. Choose = 1 for all k = l,...,n. For a set of interpolation points 
w±, W2, ■ ■ ■ , w n in the open unit disk, a set of direction vectors u\, U2, . . . , u n £ C p of length one, 
and a coordinate chart (W, i?) for the set U p of unitary p x p matrices, we define a chart (V, ip) 
by its domain 

V( Wl , w 2 , ...,w n , Ul ,u 2 ,..., u n ,W) ={Ge£ p n \ \\G^(l/w k )u k \\ < 1, G<°> G W}, 
and its coordinate map 

G -> (^i,t;2,...,t;„,^(G( ))) 

where, for k = l,...,n, the lossless functions G^ k > and the Schur parameter vectors v k are 
recursively defined in Theorem 15.51 in which we choose 

Qk = @(uk,v k ,Wk). 



Theorem 7.1 The family (V, ip) defines a C^-atlas on Cn. 
The proof is left to the reader; see also [l] and [10] . 

In approximation problems, the Douglas-Shapiro-Shields factorization is often used to represent 
a transfer function. Since in such problems the lossless factor is only determined up to a unitary 
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constant matrix, we are interested in a parametrization of the quotient space of lossless functions 
by unitary matrices. To be precise, two lossless functions G\{z) and G^z) in this quotient are 
equivalent if there exists a unitary matrix X such that G\{z) = XG 2 {z). To deal with such an 
equivalence relation, we have the identity 

which is satisfied for every lossless G and every unitary X. Therefore, if G(z) has a set of Schur 
parameters (t>i, . . . ,v n ) and a unitary matrix G^ = Dq in a chart defined by given sequences 
(wi, . . . ,w n ) of interpolation points and (ui, . . . ,u n ) of direction vectors, then X G(z) has the 
Schur parameters (X vi, . . . ,Xv n ) and the unitary matrix G^ = X Dq in the same chart. An 
atlas of the quotient space can therefore be obtained by fixing the unitary matrix G^ in each 
chart instead of letting it vary. For wi, w 2 , . . . , w n in the open unit disk, u\, u 2 , . . . , u n G C p of 
unit length and Dq, a fixed p x p unitary matrix, we define the chart (V, <p) by its domain 

V(wi,w 2 , ■ ■ ■ ,w n ,u 1 ,u 2 , ■ ■ ■ ,u n , Dq) = {G e CF n \ (l/w k )u k \\ < 1, 

and its coordinate map 

(vi,V2,...,v n ) 

where the lossless functions G^ and the Schur parameter vectors v k are again those of Theorem 
E51 in which <d k = Q(uk,v k) w k )- 

Theorem 7.2 The family (V,</5) defines a C 00 '-atlas on the quotient space C%ijlA v . 

Moreover, a unique balanced realization computed from the Schur parameters corresponds to 
every lossless function in the domain V of a chart, by iterating formula ()34|) in which U and V 
are the unitary matrices given in Theorem 16.41 These 'Schur balanced realizations' are really 
in a (local) canonical form. With each (A, B,C, D) this method associates (for a given choice 
of interpolation points and direction vectors) a unique realization (A,B,C,D). Attaching a 
canonical form with each chart, we obtain a set of overlapping canonical forms. This generalizes 
the results of [14] to the multivariable case and opens up possibilities for multivariable stable 
all-pass model reduction and approximation methods. 

Note that, except for particular choices of the interpolation points and direction vectors, these 
overlapping canonical forms in general cannot be characterized in a simple way, and cannot be 
computed directly from a given realization. However, if the interpolation points w k are the poles 
of G(z) with the direction vectors u k spanning the associated kernels of the singular matrices 
G(l/w k ), then the Schur vectors v k are all zero and the Schur balanced realization is in Schur 
form (i.e., the matrix A is triangular). If the interpolation points w k are all equal to zero and 
the vectors u k are all equal to e\ = [1, 0, . . . , 0]*, then the Schur balanced realization is positive 
upper Hessenberg. 

The overlapping balanced canonical forms for stable all-pass systems obtained in the way de- 
scribed above, also give rise to output-normal canonical forms, resp. input-normal canonical 
forms, as follows. Consider an arbitrary (asymptotically) stable system in minimal state-space 



29 



form, given by (A, B,C, D). It is not hard to show that one can find a pair (B,D) such that 
(A,B,C,D) is stable all-pass. One can then choose a state-space transformation that brings 
this system into one of the balanced canonical forms presented in this paper, which we de- 
note by (Ab, Bb,Cb, Db). The resulting state-space transformation T and the resulting pair 
(Ab,Cb) = (TAT~ 1 ,CT~ 1 ) are independent of the choice of B,D. 

Therefore applying the same state-space transformation to (A, B, C, D) one obtains a canonical 
form (A n ,B n ,C n ,D n ), where A n = A b , C n = C b , B n = TB, D n = D. Obviously the ob- 
servability Gramian of this canonical form is the identity matrix, hence it is an output-normal 
canonical form. By duality one can also construct input-normal canonical forms in the same 
fashion. These output-normal and input-normal canonical forms are of importance in model 
order reduction and system identification problems. An analogous canonical form in continuous 
time was used in [13] to construct an algorithm to find the stable system of order k, k < n, 
which is closest in the i^-norm to a given stable system of order n. In system identification 
output-normal canonical forms are used in the context of the so-called separable least squares 
approach (cf., e.g., [6]). 

The parametrizations described in this paper have been implemented in a software program 
named RARL2, dedicated to the rational approximation of multivariable discrete-time transfer 
functions in I? norm. Just as the software program Hyperion, it is based on the algorithm 
described in [10J. However, in contrast to Hyperion it has the particularity to deal with state- 
space formulations. This state-space approach is of great interest due to the good numerical 
behavior of the recursive construction of balanced realizations presented in this paper. Presently 
the main domain of application of RARL2 is the synthesis and identification of filters from partial 
frequency data. 
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